
void getChebyIndex_Anisotrop_size(int npt0, int ndim, int **icheby_iso, int vector_levels[], int vector_good[], int &npt_aniso){
	
	int points_dimensions[ndim]; 
	for (int id = 0; id<ndim; id++){
		
   		int aux = vector_levels[id]; //... approximation level for each dimension
     	
		if (aux == 0){                      // % If the approximation level in
                                            // % the i-th dimension is 0, ...
        points_dimensions[id] = 1;      	// % The number of unidimensional 
                                            // % elements is 1
     	} else {                            // % If the approximation level in
                                            // % the i-th dimension is not 0,...
        points_dimensions[id] = pow(2.0, aux)+1;  // % Compute the number of unidimensional
                                            // % elements using the formula
    	}
    	
//    	printf("\n id %d, points %3d ", id, points_dimensions[id]); 
    }
	
	for (int ip=0; ip<npt0; ip++){

//    	printf("\n ip %d \n", ip); 

		vector_good[ip] = 1; 
		for (int idim =0; idim<ndim; idim++){
			
//			printf("\t %d : %d ", icheby_iso[ip][idim], points_dimensions[idim]); 
			
			if ( icheby_iso[ip][idim] > points_dimensions[idim] ) vector_good[ip] = 0; 
		
		}
		
//		printf("\t %d ", vector_good[ip]); 
	}//ip
		
	npt_aniso = my_sum(npt0, vector_good);

}


int **getChebyIndex_Anisotrop(int npt0, int ndim, int **icheby_iso, int npt_aniso, int vector_good[]){

	int **icheby_aniso; 
	icheby_aniso = new int*[npt_aniso]; 
	for(int i = 0; i < npt_aniso; i++) icheby_aniso[i] = new int [ndim]; 
	
		int im_aniso=0; 
		for (int ip=0; ip<npt0; ip++){
			
			if (vector_good[ip] == 1){ 
				for (int idim =0; idim<ndim; idim++) icheby_aniso[im_aniso][idim] = icheby_iso[ip][idim];
				im_aniso++; 
			}
		
		}//ip	
	
	if (im_aniso != npt_aniso){
		printf("\n Errors in getChebyIndex_Anisotrop(), im_aniso = %d, npt_aniso = %d, npt0 = %d \n", im_aniso, npt_aniso, npt0); 	
	}

	return icheby_aniso; 

}



int Smolyak_Grid(const int d, const int mu, const int npt, int **chebyindex, double grid_point[]){
// output: grid_point
// input:  chebyindex = smol_elem + 1 
	
// npt is the number of rows in smol_elem

// 1. Compute the vector of extrema of Chebyshev polynomials corresponding 
// to the given level of Smolyak approximation mu
// -----------------------------------------------------------------------

// These points will be ordered as in Section 2.2.1 of JMMV(2014); e.g., for
// mu=1, the set of points is {0,-1,1}

 double pi = 3.141592653589793;
  
 int num_points_1d = pow(2, mu) + 1; 	
         
 double points_1d [num_points_1d];

 int i_max = mu+1;                 // The maximum subindex of unidimensional
                                   // set A_i whose points are used to
                                   // construct Smolyak grid of the given mu; 
                                   //  e.g., for mu=1, we consider up to 
                                   // A_i_max={-1,1} where i_max=1+1=2
				                   // A subindex of a unidimensional set of
                                   // points     
 

 int icount = 0; 
 for (int i = 1; i <= i_max; i++) {
                                                               
     // Compute the number of elements, m(i),(using m(i)=2^(i-1)+1) in the  
     // i-th unidimensional set of points; see Section 2.2.1 in JMMV (2014)
     // ---------------------------------------------------------------------
      int m_i; 
     if (i==1){
     	m_i = 1;          				// If i=1, then m(i)=1
     } else {
     	m_i = pow(2, i-1) + 1; 		// If i>1, then m(i) = 2^(i-1)+1
     }                              
     
     // Construct the extrema of Chebyshev polynomials used as unidimensional 
     // grid points in the Smolyak method
     // ---------------------------------------------------------------------
     
     double extrem_Cheb_1d; 
     
     if (m_i==1){
        
        extrem_Cheb_1d = 0;
        
        points_1d[0] = extrem_Cheb_1d; 	icount++; 
        
     } else { 
		
		for (int j = 1; j <= m_i; j++){       
		
        extrem_Cheb_1d = -cos(pi*(j-1)/(m_i-1));    
                                   // Chebyshev polynomials are defined in 
                                   // the interval [-1,1]
                                   
        if ( fabs(extrem_Cheb_1d) < 1.0e-12 ) extrem_Cheb_1d = 0.0; 
 
                                   // Round "extrem_Cheb_1d" to 0 if its  
                                   // absolute value is smaller than 1d-12
                                   
        if ( 1.0-extrem_Cheb_1d < 1.0e-12 ) extrem_Cheb_1d = 1.0; 
     
                                   // Round "extrem_Cheb_1d" to 1 if    
                                   // 1-extrem_Cheb_1d is smaller than 1d-12
                                   
        if ( 1+extrem_Cheb_1d < 1.0e-12 ) extrem_Cheb_1d = -1.0; 
    
                                   // Round "extrem_Cheb_1d" to -1 if   
                                   // 1+extrem_Cheb_1d is smaller than 1d-12
        
 //       printf("\n icount %5d, j %3d, extrem_cheb %f ", icount, j, extrem_Cheb_1d); 
		        	 	
        	 	// Check if the grid points is new
				int new_point = 1; 

				for (int i_prev =0; i_prev < icount; i_prev++){ 

					double temp = fabs( extrem_Cheb_1d - points_1d[i_prev] );
										
//			        printf("\n\t icount %5d, j %3d, extrem_cheb %f, i_prev %d, points_1d %f \t temp %f ", icount, j, extrem_Cheb_1d, i_prev, points_1d[i_prev], temp); 

					if (temp < 1.0e-10) {
						new_point = 0; // false
						break; 	
					}
			
				}
				
				if ( new_point == 1 ){
					points_1d[icount] = extrem_Cheb_1d; 
					icount++;
				}
        
  //      printf("\n icount %5d, j %3d, extrem_cheb %f ", icount, j, extrem_Cheb_1d); 
        
        
        } // j                                
     
     }
	 


} // i    

	#ifdef DEBUG_smolyak 
	for (int i_prev =0; i_prev < icount; i_prev++) printf("\n\t i_prev %d, %f ", i_prev, points_1d[i_prev]); 
	#endif 		

	if (icount != num_points_1d) { printf("\n errors in Smolyak_Grid: icount %3d, num_points_1d %3d \n", icount, num_points_1d); exit(1); }           

 
 // 2. Construct the matrix multidimensional Smolyak grid points for the   
 // required level of Smolyak approximation, mu; see JMMV (2014), Sections 2.2.3 
 // for examples
 // -------------------------------------------------------------------------
 double smol_grid[npt][d];
 for (int ip = 0; ip < npt; ip++){     // For each multidimensional grid point, ...    
     
     	#ifdef DEBUG_smolyak
		printf("\n ipt %4d, ipoly %2d \t\t ", ip, my_sum(d, chebyindex[ip]));
		#endif 
		
     for (int idim =0; idim<d; idim++){ 
     
     	int n = chebyindex[ip][idim]-1; 
     
                               // For each dimension (state variable)
                               // n denotes the subindex of the unidimensional
                               // grid point ip
        
        smol_grid[ip][idim] = points_1d[n];
                               // Find the corresponding unidimensional grid
                               // point in the vector "points_1d"
     	     	
     	grid_point[idim + ip * d] = smol_grid[ip][idim]; 
		
     	#ifdef DEBUG_smolyak
		printf(" %2d (%.4f), \t", chebyindex[ip][idim], smol_grid[ip][idim]);
		#endif
	     
     }
	      
 } // ip
	    
	    #ifdef DEBUG_smolyak
		printf("\n");
		#endif

	
	return 0; 
}



int **get_initCheby(int nrow, int ndim, int m[]){
// ilevel Chebyshev polynomials		
	int icoef=0; int **temp_icheby; 
	
	temp_icheby = new int*[nrow]; 
	for(int i = 0; i < nrow; i++) temp_icheby[i] = new int [ndim]; 
		
	if (ndim > 8) std::cerr << "\n get_initCheby(): ndim = " << ndim << " > 8 \n\n";
	
	switch(ndim){
		case 8: 
		
			for (int j0=1; j0<=m[0]; j0++){
			for (int j1=1; j1<=m[1]; j1++){
			for (int j2=1; j2<=m[2]; j2++){
			for (int j3=1; j3<=m[3]; j3++){
			
			for (int j4=1; j4<=m[4]; j4++){
			for (int j5=1; j5<=m[5]; j5++){
			for (int j6=1; j6<=m[6]; j6++){
			for (int j7=1; j7<=m[7]; j7++){
				
				temp_icheby[icoef][0] = j0;
				temp_icheby[icoef][1] = j1;
				temp_icheby[icoef][2] = j2;
				temp_icheby[icoef][3] = j3;
				
				temp_icheby[icoef][4] = j4;
				temp_icheby[icoef][5] = j5;
				temp_icheby[icoef][6] = j6;
				temp_icheby[icoef][7] = j7;
							
				icoef++; 
			
			}}}}	
			}}}}
			
			break;

		case 7: 
		
			for (int j0=1; j0<=m[0]; j0++){
			for (int j1=1; j1<=m[1]; j1++){
			for (int j2=1; j2<=m[2]; j2++){
			for (int j3=1; j3<=m[3]; j3++){
			
			for (int j4=1; j4<=m[4]; j4++){
			for (int j5=1; j5<=m[5]; j5++){
			for (int j6=1; j6<=m[6]; j6++){
				
				temp_icheby[icoef][0] = j0;
				temp_icheby[icoef][1] = j1;
				temp_icheby[icoef][2] = j2;
				temp_icheby[icoef][3] = j3;
				
				temp_icheby[icoef][4] = j4;
				temp_icheby[icoef][5] = j5;
				temp_icheby[icoef][6] = j6;
							
				icoef++; 
			
			}}}	
			}}}}
			
			break;

		case 6: 
		
			for (int j0=1; j0<=m[0]; j0++){
			for (int j1=1; j1<=m[1]; j1++){
			for (int j2=1; j2<=m[2]; j2++){
			for (int j3=1; j3<=m[3]; j3++){
			
			for (int j4=1; j4<=m[4]; j4++){
			for (int j5=1; j5<=m[5]; j5++){
				
				temp_icheby[icoef][0] = j0;
				temp_icheby[icoef][1] = j1;
				temp_icheby[icoef][2] = j2;
				temp_icheby[icoef][3] = j3;
				
				temp_icheby[icoef][4] = j4;
				temp_icheby[icoef][5] = j5;
							
				icoef++; 
			
			}}	
			}}}}
			
			break;	
	
		case 5: 
		
			for (int j0=1; j0<=m[0]; j0++){
			for (int j1=1; j1<=m[1]; j1++){
			for (int j2=1; j2<=m[2]; j2++){
			for (int j3=1; j3<=m[3]; j3++){
			
			for (int j4=1; j4<=m[4]; j4++){
				
				temp_icheby[icoef][0] = j0;
				temp_icheby[icoef][1] = j1;
				temp_icheby[icoef][2] = j2;
				temp_icheby[icoef][3] = j3;
				
				temp_icheby[icoef][4] = j4;
							
				icoef++; 
			
			}	
			}}}}
			
			break;	
	
		case 4: 
		
			for (int j0=1; j0<=m[0]; j0++){
			for (int j1=1; j1<=m[1]; j1++){
			for (int j2=1; j2<=m[2]; j2++){
			for (int j3=1; j3<=m[3]; j3++){
				temp_icheby[icoef][0] = j0;
				temp_icheby[icoef][1] = j1;
				temp_icheby[icoef][2] = j2;
				temp_icheby[icoef][3] = j3;
				icoef++; 
			}}}}	
		
			break;
			
		case 3: 
			for (int j0=1; j0<=m[0]; j0++){
			for (int j1=1; j1<=m[1]; j1++){
			for (int j2=1; j2<=m[2]; j2++){
				temp_icheby[icoef][0] = j0;
				temp_icheby[icoef][1] = j1;
				temp_icheby[icoef][2] = j2;
				icoef++; 
			}}}
		
			break;
		
		case 2: 
			for (int j0=1; j0<=m[0]; j0++){
			for (int j1=1; j1<=m[1]; j1++){
				temp_icheby[icoef][0] = j0;
				temp_icheby[icoef][1] = j1;
				icoef++; 
			}}
		
			break;

		case 1: 
			for (int j0=1; j0<=m[0]; j0++){
				temp_icheby[icoef][0] = j0;
				icoef++; 
			}
		
			break;			
		
		default: std::cerr << "Unknown Dimensional Type";
	}
	

	return temp_icheby; 
}



int **getChebyIndex(int l_min, int l_max, int npt, int ndim){
// npt: the number of distinct set of tensor-product multivariate chebyshev polynomial

// P^I = p^(i_1) p^(i_2) ... p^(i_d), and i_1 + i_2 + ... + i_d = I \in [ level_min, level_max ]
//   where p^i is a univariate polynomials of degree m_i = 2^(i-1)+1

// The combination of tensor-product multivariate polynomial
//     T_(l_1) * T_(l_2) * ... * T_(l_dim), ChebyIndex[.][l_d] = l_d
// The number of such unique tensor-product multivariate polynomial is npt. 
// Therefore, out_chebyindex has npt rows and dim column. 
	
	int **icheby; 
	icheby = new int*[npt]; 
	for(int i = 0; i < npt; i++) icheby[i] = new int [ndim]; 
	
	
	int *level_1d = new int [ndim];
	bool more = false;
    int h = 0;
    int t = 0;
    
    for (int id = 0; id<ndim; id++)  icheby[0][id] = 0; 
	int im = 0; 
	for (int l = l_min; l <= l_max; l ++){
	
    	for ( ; ; )  {
        // compute the compositions of the integer N=level into K=dim_num parts. 
      
        comp_next ( l, ndim, level_1d, &more, &h, &t );
					 
			int m[ndim]; int n_m=1; 
			for (int id =0; id<ndim; id++){
				if (level_1d[id]==0){	// id=level_1d[id]+1=1
					m[id] = 1; 
				} else {
					m[id] = (int) pow(2.0, level_1d[id]) + 1;
				}
				n_m *= m[id];	// multivariate polynomials
			}
			
			int **temp_icheby; 
			temp_icheby = new int* [n_m]; 
			for(int i = 0; i < n_m; i++) temp_icheby[i] = new int [ndim]; 
			temp_icheby = get_initCheby(n_m, ndim, m);
		
					
			// Check if the combination of Chebyshev polynomials is new
			for (int in=0; in<n_m; in++){
				int new_poly = 1; 
				for (int i_prev =0; i_prev < im; i_prev++){
					int temp = 0; 
					for (int id = 0; id<ndim; id++){
						temp += abs( temp_icheby[in][id] - icheby[i_prev][id] );

					}
					if (temp == 0) {
						new_poly = 0; // false
						break; 	
					}				
				}
				
				// Update the combination of chebyshev polynomial into ChebyIndex if new
				if (new_poly==1) {

					for (int id = 0; id<ndim; id++) {
					
						icheby[im][id] = temp_icheby[in][id]; 
					
					}
					im++;
				}

			}
	
			for(int i = 0; i < n_m; i++) delete [] temp_icheby[i];
				
		if ( !more ){
        	break;
        }
		
	} 

	} // l 
	
	
	if (im != npt){
		printf("\n Errors in getChebyIndex(), im = %d, npt = %d \n", im, npt); 	
	}

	return icheby; 

}




struct struct_smolyak{

	int order, npt, ndim; 

	struct_smolyak() : order(ORDER_smolyak), npt(NUM_coefEmax), ndim(DIM_smolyak) { 
		
	}
		
	void ChebyPoly(int **chebyind, double state[], gsl_vector * tensorCheby); 	

	int  ChebyReg(int **chebyind, double *y, double **X_states, int N, gsl_vector *ols_coef); 
	
	double Interpolate(int **chebyind, double state[], gsl_vector * coef); 
		
	void PolyInvert (int **chebyind, double **X_states, gsl_matrix * mat_polyinv); 
	void PolyCoef(int **chebyind, gsl_matrix *polyinv, double *y, gsl_vector *coef); 
	
	double PolyInterp(int **chebyind, double state[], double state_min[], double state_max[], double coef[]); 

	double Derivative(int **chebyind, double state[], double state_min[], double state_max[], double coef[], const int ideriv); 
	
}; 


void   struct_smolyak::ChebyPoly(int **chebyind, double state[], gsl_vector * tensorCheby){
	double tempcheby[order][ndim]; 

	for (int id =0; id< ndim; id++){
		tempcheby[0][id] = 1.0; 
		tempcheby[1][id] = state[id];
		
		for(int j=2; j < order ; j++)
			tempcheby[j][id] = 2.0*state[id]*tempcheby[j-1][id] - tempcheby[j-2][id];
	}

	gsl_vector_set_all(tensorCheby, 1.0);
	
	for (int i = 0; i< npt; i++){

		double temp_t = 1.0;
		for (int id=0; id<ndim; id++){
			int icheby = 0; 
			icheby = chebyind[i][id] - 1; 
			temp_t *= tempcheby[icheby][id]; 
		}
		
		gsl_vector_set(tensorCheby, i, temp_t); 
	}
	
}


int    struct_smolyak::ChebyReg(int **chebyind, double *y, double **states_X, int N, gsl_vector *ols_coef){
	gsl_vector *ols_y		= gsl_vector_alloc(N); 
	gsl_vector *ols_tensor 	= gsl_vector_alloc(N); 
	gsl_matrix *ols_x 		= gsl_matrix_alloc(N, N);	
	gsl_matrix *ols_cov  	= gsl_matrix_alloc(N, N); 
	
	
	for (int ipt = 0; ipt< npt; ipt++){
		
		gsl_vector_set(ols_y, ipt, y[ipt]); 
				
		ChebyPoly(chebyind, states_X[ipt], ols_tensor); 
		
		gsl_matrix_set_row(ols_x, ipt, ols_tensor);
	}
	gsl_multifit_linear_workspace * work =  gsl_multifit_linear_alloc (N, N);	
	double chisq;	
	
	gsl_vector_set_zero(ols_coef); 
	
    gsl_multifit_linear (ols_x, ols_y, ols_coef, ols_cov, &chisq, work);
	
    gsl_multifit_linear_free (work);
    gsl_vector_free(ols_tensor); 
	gsl_matrix_free(ols_x); gsl_vector_free(ols_y); gsl_matrix_free(ols_cov);
    
    return 0 ; 
	
}


double struct_smolyak::Interpolate(int **chebyind, double state[], gsl_vector * coef){
	gsl_vector * tensorCheby = gsl_vector_alloc(npt);	
	
	ChebyPoly(chebyind, state, tensorCheby); 	
	
	double v=0.0; 
	for (int ipt = 0; ipt< npt; ipt++)
		v += gsl_vector_get(coef, ipt)*gsl_vector_get(tensorCheby, ipt) ; 
	
	gsl_vector_free(tensorCheby);	
	
	return v; 

}


void   struct_smolyak::PolyInvert(int **chebyind, double **X_states, gsl_matrix * mat_polyinv){

	gsl_matrix * matX        = gsl_matrix_alloc(npt, npt); 
	gsl_vector * tensorCheby = gsl_vector_alloc(npt);
	gsl_matrix * matXX       = gsl_matrix_alloc(npt, npt); 
	gsl_matrix * invXX       = gsl_matrix_alloc(npt, npt); 

	for (int ipt = 0; ipt < npt; ipt++){

		ChebyPoly(chebyind, X_states[ipt], tensorCheby); 				
		
		gsl_matrix_set_row( matX, ipt, tensorCheby);

	}
	
    gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, matX, matX, 0.0, matXX);
	
	int s;	
	gsl_permutation * p = gsl_permutation_alloc (npt);	
	
	gsl_linalg_LU_decomp (matXX, p, &s); 
	
	gsl_linalg_LU_invert (matXX, p, invXX); 
	
    gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1.0, invXX, matX, 0.0, mat_polyinv);
	
	gsl_matrix_free(matX); 
	gsl_vector_free(tensorCheby); 	
	gsl_matrix_free(matXX); 
	gsl_matrix_free(invXX); 
	
}


void   struct_smolyak::PolyCoef(int **chebyind, gsl_matrix *polyinv, double *y, gsl_vector *coef){
	gsl_vector_set_zero(coef); 
	
	for (int i1 = 0; i1 < npt; i1 ++){
		
		double s = 0.0; 
		for (int i2=0; i2<npt; i2++){
			s += gsl_matrix_get(polyinv, i1, i2) * y[i2]; 
		}	
		gsl_vector_set(coef, i1, s); 
		
	}
	
}	


double struct_smolyak::PolyInterp(int **chebyind, double state[], double state_min[], double state_max[], double coef[]){

	#ifdef DEBUG 
		
		for (int idim=0; idim<ndim; idim++){
			if ( isnan(state[idim]) || state[idim] < state_min[idim] - 1.0e-10 || state[idim] > state_max[idim] + 1.0e-10 ){
				printf("\n\n Errors in struct_smolyak::PolyInterp: idim %2d, state %f, [%f, %f] \n ... exiting ... \n", idim, state[idim], state_min[idim], state_max[idim] ); 
				exit(1); 
			}
		}	
		
	#endif 
	
	
	double state_norm[ndim]; 
	
	for (int id=0; id<ndim; id++)
		state_norm[id] = 2.0*(state[id]-state_min[id])/(state_max[id]-state_min[id]) - 1.0;
			
	gsl_vector * tensorCheby = gsl_vector_alloc(npt);	
	
	ChebyPoly(chebyind, state_norm, tensorCheby); 	
	
	double v=0.0; 
	for (int ipt = 0; ipt< npt; ipt++)
		v += coef[ipt]*gsl_vector_get(tensorCheby, ipt) ; 
	
	gsl_vector_free(tensorCheby);	
	
	return v; 
		
}



double struct_smolyak::Derivative(int **chebyind, double state1[], double state_min[], double state_max[], double coef[], const int iderive){

	#ifdef DEBUG 
		
		for (int idim=0; idim<ndim; idim++){
			if ( isnan(state1[idim]) || state1[idim] < state_min[idim] - 1.0e-10 || state1[idim] > state_max[idim] + 1.0e-10 ){
				printf("\n\n Errors in struct_smolyak::PolyInterp: idim %2d, state %f, [%f, %f] \n ... exiting ... \n", idim, state1[idim], state_min[idim], state_max[idim] ); 
				exit(1); 
			}
		}	
		
	#endif 
	
	
	double state_norm[ndim]; 
	
 	for (int id=0; id<ndim; id++)
		state_norm[id] = 2.0*(state1[id]-state_min[id])/(state_max[id]-state_min[id]) - 1.0;
	
	
	gsl_vector * tensorCheby = gsl_vector_alloc(npt);	
	
	
	double tempcheby[order][ndim];  
	double cheby2[order][ndim];  
	
 	for (int id =0; id< ndim; id++){
		
 		tempcheby[0][id] = 1.0; 
		tempcheby[1][id] = state_norm[id];
		
		for(int j=2; j < order ; j++)
			tempcheby[j][id] = 2.0*state_norm[id]*tempcheby[j-1][id] - tempcheby[j-2][id];
		
		
		// ... The Chebyshev polynomials of the second kind U_n, d(T_n)/dx = n * U_{n-1}
 		cheby2[0][id] = 1.0; 
		cheby2[1][id] = 2.0*state_norm[id];
		
		for(int j=2; j < order ; j++)
			cheby2[j][id] = 2.0*state_norm[id]*cheby2[j-1][id] - cheby2[j-2][id];
			
	}
		
	gsl_vector_set_all(tensorCheby, 1.0);
	
 	for (int i = 0; i< npt; i++){

		double temp_t = 1.0;
		for (int id=0; id<ndim; id++){
			int icheby = 0; 
			icheby = chebyind[i][id] - 1;  
			
			if (id == iderive ){ 
			 	if (icheby >=1){ 
				temp_t *= (icheby * cheby2[icheby-1][id]); 
				} else { temp_t = 0.0; }
			} else {
				temp_t *= tempcheby[icheby][id]; 
			}
		
		}
			
		gsl_vector_set(tensorCheby, i, temp_t); 
	}
	
		

	double v=0.0; 
	for (int ipt = 0; ipt< npt; ipt++)
		v += coef[ipt]*gsl_vector_get(tensorCheby, ipt) ; 
	
	gsl_vector_free(tensorCheby);	
	
	return (v * 2.0/(state_max[iderive]-state_min[iderive])); 
		
}





double ChebyReg_testfn(double temp_state[]){
	
//	return ( 1.0/ ( 1.0 - 0.99*exp( 0.002*temp_state[0] + 0.005*temp_state[1])) * (pow( temp_state[2]+1.0 +1.0, 1.0-cu_g ) - 1.0 ) / (1.0 - cu_g) ); 

	return ( 3.0 * temp_state[2] + 4.0 *temp_state[2]*temp_state[2] + 8.0 *temp_state[2]*temp_state[2]*temp_state[2] + exp(temp_state[2])); 

}

double ChebyReg_testfn_der(double temp_state[]){

	double	der = 0.0; 
	
//		der = 1.0/ ( 1.0 - 0.99*exp( 0.002*temp_state[0] + 0.005*temp_state[1])) * pow( temp_state[2]+1.0 +1.0, -cu_g ); 

	der = 3.0 + 8.0 * temp_state[2] + 24.0*temp_state[2]*temp_state[2] + exp(temp_state[2]); 
	
	return der; 	
			
}


int ChebyReg_norm_test(double **states, double **states_norm, double state_min[], double state_max[]){
	
	time_t start_t,end_t; double diff_t; 
	time (&start_t);
	
	gsl_vector * coefValue = gsl_vector_alloc(NUM_coefEmax); 
	double value[NUM_coefEmax]; double der[NUM_coefEmax]; 
	
	for (int ipt =0 ; ipt< NUM_coefEmax; ipt++){
		value[ipt] = ChebyReg_testfn(states[ipt]);
		der[ipt] = ChebyReg_testfn_der(states[ipt]);
	}
	
	
	struct_smolyak smolyak; 
		
		// Update coefficient
//		smolyak.PolyInvert(states_norm, POLYINV_smolyak); 	
		
		smolyak.PolyCoef(ChebyIndex, POLYINV_smolyak, value, coefValue); 	

		
		double coef[NUM_coefEmax]; 
		for (int ipt =0 ; ipt< NUM_coefEmax; ipt++) coef[ipt] = gsl_vector_get(coefValue, ipt);		
		
		
		// Interpolation
		
		printf( "\n\n TEST ON THE GRID \n\n"); 
		
		double abserror = 0.0; 			
		for (int ipt =0 ; ipt< NUM_coefEmax; ipt++)
		{
			double vhat = smolyak.PolyInterp(ChebyIndex, states[ipt], state_min, state_max, coef); 
			double derhat = smolyak.Derivative(ChebyIndex, states[ipt], state_min, state_max, coef, 2); 
			printf("\n ipt %3d, value %f, \t vhat %f \t s %f \t der %f, \t derhat %f", ipt, value[ipt], vhat, states[ipt][2], der[ipt], derhat); 
			abserror += fabs ( vhat - value[ipt] ); 
		} 
			printf( "\n \t mean(emax) = %f, \t abserror = %f \n", my_mean(NUM_coefEmax, value), abserror/NUM_coefEmax); 
						
		double temp_state[] = {states[0][0], states[0][1], 0.2, states[0][3], states[0][4], states[0][5]}; // {0.5, -0.5, 0.2, 0.5, 1.0};	
		
		double temp_value = ChebyReg_testfn(temp_state); 
				
		double vhat = smolyak.PolyInterp(ChebyIndex, temp_state, state_min, state_max, coef); 
			
		printf( "\n\n haha: true = %f, \t interp = %f \n", temp_value, vhat ); 
	
	
	time (&end_t);	diff_t = difftime (end_t, start_t);
	printf ("\n\n It took you %.10f seconds to interpolation \n", diff_t );		
	
	return 0; 
}

/*
int ChebyReg_test(double **states){

	time_t start_t,end_t; double diff_t; 
	time (&start_t);
	
	gsl_vector * coefValue = gsl_vector_alloc(NUM_coefEmax); 
	double value[NUM_coefEmax]; 	
	
	for (int ipt =0 ; ipt< NUM_coefEmax; ipt++){
		value[ipt] = ChebyReg_testfn(states[ipt]);
	}
	
	struct_smolyak smolyak; 
		
		// Update coefficient
		smolyak.ChebyReg(value, states, NUM_coefEmax, coefValue);		
		
		// Interpolation
		
		printf( "\n\n TEST ON THE GRID \n\n"); 	
		double abserror = 0.0; 		
		for (int ipt =0 ; ipt< NUM_coefEmax; ipt++)
		{	
			double vhat = smolyak.Interpolate(states[ipt], coefValue); 
			
//			printf( "\n\n ipt %d, true = %f, \t interp = %f, \t state: ", ipt, value[ipt], vhat); 
//			for (int is=0; is<DIM_smolyak; is++) printf(" %f", states[ipt][is]); 
			
			abserror += fabs ( vhat - value[ipt] ); 
		} 
		printf( " \t mean(emax) = %f, \t abserror = %f \n", my_mean(NUM_coefEmax, value), abserror/NUM_coefEmax); 
		
		double temp_state[] = {0.5, -0.5, 1.0, 100.0, 0.5, 0.0, 12.0}; 
		double temp_value = ChebyReg_testfn(temp_state); 
		
		double vhat = smolyak.Interpolate(temp_state, coefValue); 
			
		printf( "\n\n haha: true = %f, \t interp = %f \n", temp_value, vhat ); 

				
	time (&end_t);	diff_t = difftime (end_t, start_t);
	printf ("\n\n It took you %.10f seconds to interpolation \n", diff_t );		
	
	
	return 0; 

}
*/

